Explaining roe deer movement in the Italian alps

This project investigates movement of GPS collard roe deer in the area of Trento, Italy. It compares recorded GPS data of five animals (Agostino, Alessandra, Daniela, Decimo, and Sandro) with simulated roe deer trajectories in order to assess to what extent roe deer movement can be explained by altitude and/or land use. The results from this project show that the roe deer stick to a relatively small home range and stable altitude compared to the simulated locations. Roe deer had a preference for forest and pasture. However, the difference with simulated data is small, which can be explained by the high prevalence of these land use types in the home range of the roe deer.

GPS-collared roe deer

Prepare working environment

# check if libraries are installed
if (!require("dplyr")) install.packages("dplyr")
if (!require("sp")) install.packages("sp")
if (!require("sf")) install.packages("sf")
if (!require("DBI")) install.packages("DBI")
if (!require("RSQLite")) install.packages("RSQLite")
if (!require("raster")) install.packages("raster")
if (!require("adehabitatLT")) install.packages("adehabitatLT")
if (!require("lubridate")) install.packages("lubridate")
if (!require("ggplot2")) install.packages("ggplot2")
if (!require("gridExtra")) install.packages("gridExtra")
if (!require("leaflet")) install.packages("leaflet")
if (!require("rlist")) install.packages("rlist")
if (!require("rasterVis")) install.packages("rasterVis")

# load libraries
library(dplyr)
library(sp)
library(sf)
library(DBI)
library(raster)
library(adehabitatLT)
library(lubridate)
library(ggplot2)
library(gridExtra)
library(leaflet)
library(rlist)
library(rasterVis)

# import functions
source("R/retrieveData.R")
source("R/saveDataInDB.R")
source("R/simulateMovement.R")
source("R/createSpatialLine.R")
source("R/extractRasterValues.R")
source("R/plotProfiles.R")
source("R/combineSimTrack.R")
source("R/sumStats.R")
source("R/boxplotProfiles.R")
source("R/addLegend.R")

# project variables
# downloadlinks
movementDataURL = "extras.springer.com/2014/978-3-319-03742-4/trackingDB_datasets.zip"
environmentalURL = "https://drive.google.com/uc?export=download&id=1nb7q8TQGAkzTcGtSMS0Nig_DAhB2lClJ"
landUseLegendURL  = "https://www.eea.europa.eu/data-and-maps/data/corine-land-cover-2006-raster-1/corine-land-cover-classes-and/clc_legend.csv/at_download/file"
# projections
proj4WGS84 <- "+proj=longlat +datum=WGS84 +no_defs"
proj4zone32 <- "+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs"
proj4Mercator <- "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"

Obtain data

Google Earth data collection

# download movement data and unzip
folder = "data"
retrieveData(movementDataURL, folder, unzip = TRUE)

# download environmental data and unzip
retrieveData(environmentalURL, folder, unzip = TRUE)

# download legend for corine land cover (landUse)
retrieveData(landUseLegendURL, folder, filename = "landUseLegend.csv")

# Load environmental data to global environment
DEM <- raster(list.files(path = "data/GoogleEarthData", pattern = glob2rx('*DEM*.tif'), full.names = TRUE))
landsat <- brick(list.files(path = "data/GoogleEarthData", pattern = glob2rx('*Landsat*.tif'), full.names = TRUE))
landUse <- raster(list.files(path = "data/GoogleEarthData", pattern = glob2rx('*LandUse*.tif'), full.names = TRUE))
slope <- raster(list.files(path = "data/GoogleEarthData", pattern = glob2rx('*Slope*.tif'), full.names = TRUE))
landUseLegend <- read.csv(list.files(path = folder, pattern = glob2rx('*.csv'), full.names = TRUE))

Database management

We use an example dataset of roe deer movement from a book about spatial database management (Urbano and Cagnacci, 2014). The data are thus designed for use in a database and set up as interlinking datasets. Direct use in R is possible, but undesirable due to a lack of structure. Therefore, we set up an SQLite database in memory to process the data.

# create database connection
movementDB <- dbConnect(RSQLite::SQLite(), "")
# get files
gpsData <- list.files(path = "data/tracking_db/data/sensors_data/", pattern = glob2rx('GSM0*.csv'), full.names = TRUE)
sensorsAnimals <- list.files(path = "data/tracking_db/data/sensors_animals/", pattern = glob2rx('*.csv'), full.names = TRUE)
sensors <- list.files(path = "data/tracking_db/data/sensors/", pattern = glob2rx('*.csv'), full.names = TRUE)
animals <- list.files(path = "data/tracking_db/data/animals/", pattern = glob2rx('*.csv'), full.names = TRUE)
# specify column names for database tables
gpsDataColumnNames <- c("gpsSensorsCode","lineNo", "utcDate", "utcTime", "lmtDate", "lmtTime", "ecefX", "ecefY","ecefZ", "latitude", "longitude", 
"height", "dop", "nav", "validated", "satsUsed", "ch01SatId", "ch01SatCnr", "ch02SatId", "ch02SatCnr", "ch03SatId", "ch03SatCnr", "ch04SatId", "ch04SatCnr", "ch05SatId", "ch05SatCnr", "ch06SatId", "ch06SatCnr", "ch07SatId", "ch07SatCnr", "ch08SatId", "ch08SatCnr", "ch09SatId", "ch09SatCnr", "ch10SatId", "ch10SatCnr", "ch11SatId", "ch11SatCnr", "ch12SatId", "ch12SatCnr", "mainVol", "buVol", "temp", "easting", "northing", "remarks")
sensorsAnimalsColumnNames <- c("animalsId","gpsSensorsId","startTime","EndTime","notes")
sensorsColumnNames <-c("gpsSensorsId", "gpsSensorsCode", "purchaseDate", "frequency", "vendor","model","sim")
animalsColumnNames <- c("animalsId", "animalsCode", "name","sex","ageClassCode", "speciesCode")
# save data in database
saveDataInDB(gpsData, movementDB, "gpsData", gpsDataColumnNames, TRUE)
saveDataInDB(sensorsAnimals, movementDB, "sensorsAnimals", sensorsAnimalsColumnNames, FALSE)
saveDataInDB(sensors, movementDB, "sensors", sensorsColumnNames, FALSE)
saveDataInDB(animals, movementDB, "animals", animalsColumnNames, FALSE)

The database can be queried with SQL. Since we want to link the GPS data with specific animals we need to join these tables. However, these tables do not have any common fields that can be used for a straightforward join. Embedding of SQL statements allows us to get to the required combination of data.

# query data from database
query <- '
SELECT
animals.name,
animals.sex,
dataWithAnimalId.utcDate, dataWithAnimalId.utcTime, 
dataWithAnimalId.latitude, dataWithAnimalId.longitude,
dataWithAnimalId.height
FROM
  (SELECT 
  dataWithSensorId.gpsSensorsCode, 
  dataWithSensorId.utcDate, dataWithSensorId.utcTime, 
  dataWithSensorId.latitude, dataWithSensorId.longitude, 
  dataWithSensorId.height,
  dataWithSensorId.gpsSensorsId, 
  sensorsAnimals.animalsId 
  FROM 
    (SELECT 
    gpsData.gpsSensorsCode, 
    gpsData.utcDate, gpsData.utcTime, 
    gpsData.latitude, gpsData.longitude, 
    gpsData.height,
    sensors.gpsSensorsId 
    FROM 
    gpsData 
    LEFT JOIN 
    sensors ON gpsData.gpsSensorsCode = sensors.gpsSensorsCode) AS dataWithSensorId 
  LEFT JOIN 
  sensorsAnimals ON dataWithSensorId.gpsSensorsId = sensorsAnimals.gpsSensorsId) AS dataWithAnimalId 
LEFT JOIN
animals ON dataWithAnimalId.animalsId = animals.animalsId
;'
selectedMovementData <- dbGetQuery(movementDB, query)

Pre-processing

At this point we have the data we need. However, all data is stored as strings. The correct classes are assigned, before we can create a geometry out of the latitude and longitude columns. We choose the WGS84 zone 32 (EPSG:32632) as coordinate system as it provides enough detail for our analyses.

# correct classes
selectedMovementData$sex <- factor(selectedMovementData$sex)
class(selectedMovementData$latitude) <- "numeric"
class(selectedMovementData$longitude) <- "numeric"
class(selectedMovementData$height) <- "numeric"
selectedMovementData$date <- dmy_hms(paste(selectedMovementData$utcDate, selectedMovementData$utcTime), tz = "UTC")
# remove redundant columns
selectedMovementData[3:4] <- list(NULL)

# remove rows with NA for location
selectedMovementDataNoNA <- selectedMovementData[!is.na(selectedMovementData$latitude) | !is.na(selectedMovementData$longitude),]

# create spatial data frame projected in WGS84 zone 32
coords <- selectedMovementDataNoNA[,c(4:3)]
WGS84movementSDF <- SpatialPointsDataFrame(coords = coords, data = selectedMovementDataNoNA, proj4string = CRS(proj4WGS84))
movementSDF <- spTransform(WGS84movementSDF, CRS(proj4zone32))

# create trajectories
movementTrajectories <- as.ltraj(xy = coordinates(movementSDF), date = movementSDF$date, id = movementSDF$name)

# and create a dataframe to group our data and add height information
dataframeTrajectories <- ld(movementTrajectories)
groups <- dataframeTrajectories %>% mutate(height = movementSDF$height) %>% group_by(id)
trajectories <- dl(groups)

We are curious where “our” roe deer are going! For visualization purposes we need to create spatial lines of the GPS recorded trajectories.

# select distinct animals
animalNames <- distinct(groups, id)[[1]]
for (animal in animalNames){
  # get location data from nth animal
  locationdatainit <- groups[(groups$id == animal),]
  
  # create spatial line
  assign(paste0("trajectinit",animal), createSpatialLine(locationdata = locationdatainit, proj4dataframe = proj4zone32, proj4line = proj4WGS84, lineId = animal))
}

An interactive leaflet map allows us to explore the GPS recorded movement of our roe deer, Agostino, Alessandra, Daniela, Decimo, and Sandro.

m <- leaflet() %>% 
  setView(lng=11.0367, lat=46.0331, zoom = 11) %>%
  # baselayers
  addTiles(group = 'Street map') %>%  
  addProviderTiles('Stamen.Terrain', group = 'Terrain') %>% 
  addProviderTiles('Esri.WorldImagery', group = 'Satellite') %>%
  # overlay groups
  addPolylines(data = trajectinitAgostino, weight = 1, opacity = 0.8, color = "red", group = 'Agostino') %>% 
  addPolylines(data = trajectinitAlessandra, weight = 1, opacity = 0.8, color = "blue", group = 'Alessandra') %>% 
  addPolylines(data = trajectinitDaniela, weight = 1, opacity = 0.8, color = "orange", group = 'Daniela') %>% 
  addPolylines(data = trajectinitDecimo, weight = 1, opacity = 0.8, color = "yellow", group = 'Decimo') %>% 
  addPolylines(data = trajectinitSandro, weight = 1, opacity = 0.8, color = "purple", group = 'Sandro') %>% 
  addLayersControl(
    baseGroups = c('Street map', 'Terrain', 'Satellite'),
    overlayGroups = c('Agostino', 'Alessandra', 'Daniela', 'Decimo', 'Sandro'),
    options = layersControlOptions(collapsed = FALSE)
  )
m 

Oops! The data clearly contain some outliers. These need to be removed in order to perform reliable analyses. But how do we define which points are outliers? We look at boxplots of different variables in order to determine a suitable method to identify outliers.

# plot
p1 <- ggplot(groups, aes(x=id, y=x, fill = id)) + geom_boxplot(alpha=0.2) + theme(legend.position = "none") + ggtitle("x-coordinates of locations")
p2 <- ggplot(groups, aes(x=id, y=y, fill = id)) + geom_boxplot(alpha=0.2) + theme(legend.position = "none") + ggtitle("y-coordinates of locations")
p3 <- ggplot(groups, aes(x=id, y=dist, fill = id)) + geom_boxplot(alpha=0.2) + theme(legend.position = "none") + ggtitle("Distance between relocations")
p4 <- ggplot(groups, aes(x=id, y=abs.angle, fill = id)) + geom_boxplot(alpha=0.2) + theme(legend.position = "none") + ggtitle("Absolute angle between relocations")
p5 <- ggplot(groups, aes(x=id, y=rel.angle, fill = id)) + geom_boxplot(alpha=0.2) + theme(legend.position = "none") + ggtitle("Relative angle between relocations")
p6 <- ggplot(groups, aes(x=id, y=dt, fill = id)) + geom_boxplot(alpha=0.2) + theme(legend.position = "none") + ggtitle("Time between relocations")
grid.arrange(p1,p2,p3,p4,p5,p6, layout_matrix = rbind(c(1,2),c(3,4),c(5,6)))

These boxplots show the outliers very clearly in the x and y-coordinates as well as the distance and time between relocations. In previous steps we removed some measurements with NA’s for latitude and longitude to be able to convert the dataframe to a spatial dataframe. As a result, the time between relocations is not 100% accurate. We will therefore use only x-coordinates, y-coordinates and distance between relations to clean our dataset. Tresholds for the x and y-coordinates are based purely on these boxplots. To determine the treshold for the distance we take a look at the distribution of steplengths.

# plot histograms of distance between location points
par(mfrow=c(3, 2))
for (i in 1:length(id(movementTrajectories))){
  hist(movementTrajectories[[i]]$dist, xlim = c(0, 5000), ylim = c(0, 1500), breaks = seq(0, 1000000, by = 100), xlab = "distance (m)", main = paste("Distance between relocations by", id(movementTrajectories)[i]))
}

Looking at the distribution of steplengths and taking into account the consideration that roe deer generally do not travel great distances, a moved distance of 1500m or higher in one relocation is considered unrealistic.

# tresholds
xMin <- 650000
xMax <- 700000
yMin <- 5092500
yMax <- 5107500
distMax <- 1500

# remove outliers based on x-coordinates and distance between relocations
groupsCleaned <- groups %>% filter(x < xMax) %>% filter(x > xMin) %>% filter(y < yMax) %>% filter(y > yMin) %>% filter(dist < distMax)

# cast back to ltraj 
movementTrajectoriesCleaned <- dl(groupsCleaned)

For visualization purposes we create Spatial Lines of the recorded trajectories without outliers.

# select distinct animals
animalNames <- distinct(groupsCleaned, id)[[1]]
for (animal in animalNames){
  # get location data from nth animal
  locationdata <- groupsCleaned[(groupsCleaned$id == animal),]
  
  # create spatial line
  assign(paste0("traject",animal), createSpatialLine(locationdata = locationdata, proj4dataframe = proj4zone32, proj4line = proj4WGS84, lineId = animal))
}

We can now visualize the cleaned movement patterns of our beloved roe deer!

m <- leaflet() %>% 
  setView(lng=11.0367, lat=46.0331, zoom = 11) %>%
  # baselayers
  addTiles(group = 'Street map') %>%  
  addProviderTiles('Stamen.Terrain', group = 'Terrain') %>% 
  addProviderTiles('Esri.WorldImagery', group = 'Satellite') %>%
  # overlay groups
  addPolylines(data = trajectAgostino, weight = 1, opacity = 0.8, color = "red", group = 'Agostino') %>% 
  addPolylines(data = trajectAlessandra, weight = 1, opacity = 0.8, color = "blue", group = 'Alessandra') %>% 
  addPolylines(data = trajectDaniela, weight = 1, opacity = 0.8, color = "orange", group = 'Daniela') %>% 
  addPolylines(data = trajectDecimo, weight = 1, opacity = 0.8, color = "yellow", group = 'Decimo') %>% 
  addPolylines(data = trajectSandro, weight = 1, opacity = 0.8, color = "purple", group = 'Sandro') %>% 
  addLayersControl(
    baseGroups = c('Street map', 'Terrain', 'Satellite'),
    overlayGroups = c('Agostino', 'Alessandra', 'Daniela', 'Decimo', 'Sandro'),
    options = layersControlOptions(collapsed = FALSE)
  )
m 

Note! The data might still contain a few wrongly recorded points (we are not sure whether Agostino really visited his aunt in Stravino for a coffee). However, we are not justified to remove these potential outliers on the basis of our previously defined cut-off point.

Boxplots allow for further assessment of the cleaned data. Looks a lot better, doesn’t it?

# boxplots
p1 <- ggplot(groupsCleaned, aes(x=id, y=x, fill = id)) + geom_boxplot(alpha=0.2) + theme(legend.position = "none") + ggtitle("x-coordinates of locations")
p2 <- ggplot(groupsCleaned, aes(x=id, y=y, fill = id)) + geom_boxplot(alpha=0.2) + theme(legend.position = "none") + ggtitle("y-coordinates of locations")
p3 <- ggplot(groupsCleaned, aes(x=id, y=dist, fill = id)) + geom_boxplot(alpha=0.2) + theme(legend.position = "none") + ggtitle("Distance between relocations")
p4 <- ggplot(groupsCleaned, aes(x=id, y=abs.angle, fill = id)) + geom_boxplot(alpha=0.2) + theme(legend.position = "none") + ggtitle("Absolute angle between relocations")
p5 <- ggplot(groupsCleaned, aes(x=id, y=rel.angle, fill = id)) + geom_boxplot(alpha=0.2) + theme(legend.position = "none") + ggtitle("Relative angle between relocations")
p6 <- ggplot(groupsCleaned, aes(x=id, y=dt, fill = id)) + geom_boxplot(alpha=0.2) + theme(legend.position = "none") + ggtitle("Time between relocations")
grid.arrange(p1,p2,p3,p4,p5,p6, layout_matrix = rbind(c(1,2),c(3,4),c(5,6)))

Data validation

The data is ready for the analyses!

Before we start our movement analysis, we will assess the reliability of our data. The altitude measurement of the GPS data needs to be validated through a comparison with the DEM. GPS data are not reliable if its mean altitude differs significantly from the DEM values at corresponding spatial points.

A boxplot allows us to have a first look at the difference in height measurement between DEM and GPS.

# calculate DEM height values at point locations for which GPS height data is available.
heightDEM <- extract(DEM, as.data.frame(groupsCleaned[,c("x","y")]), method='simple', na.rm=TRUE, df=TRUE)
heightDEM$animal <- unlist(groupsCleaned[,"id"])
colnames(heightDEM) <- c("ID", "rasterVal", "animal")
heightData <- data.frame(group = rep(c("DEM", "GPS"), each = nrow(heightDEM)), height = c(heightDEM$rasterVal, groupsCleaned$height))
# visualize
boxplot(height ~ group, heightData, main = "Height boxplot GPS and DEM", ylab = "Height")

Ooh no! The GPS data contains large outliers (the alps are not even 20 000m high!). Did Allessandro also have a coffee in outer space?

A paired t-test is performed to assess the reliability the measurements.

# perform a paired t-test to compare mean GPS and DEM height
tTest <- t.test(heightDEM$rasterVal, groupsCleaned$height, paired = TRUE, alternative = "two.sided")
print("Results of the paired t-test:")

[1] "Results of the paired t-test:"

tTest

    Paired t-test

data:  heightDEM$rasterVal and groupsCleaned$height
t = -9.7866, df = 10973, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -49.59823 -33.04543
sample estimates:
mean of the differences 
              -41.32183 

The paired t-test confirms that there is a significant difference (p-value < 2.2e-16) in mean altitude between the GPS height measurement and the DEM height values. The mean difference between the two measurements amounts to 41 meters. This confirms our suspicion that the GPS altitude measurements are not always reliable. Therefore, we will use DEM altitude measurement for our further analysis.

Now that we have obtained reliable altitude measurements, we will extract landuse values for the recorded trajectories.

# extract landuse values of recorded trajectories
landuseTracks <- extract(landUse, as.data.frame(groupsCleaned[,c("x","y")]), method='simple', na.rm=TRUE, df=TRUE)
landuseTracks$animal <- unlist(groupsCleaned[,"id"])
colnames(landuseTracks) <- c("ID", "rasterVal", "animal")

Simulations

In the next section we will start with the simulations of animal trajectories. We will use the simm.crw function of the adehabitatLT package for the simulations. We wrapped this function in the simulateMovement function so that it is easy to specify the number of simulations and save all simulation in an convenient dataformat. The simulations start at the same location as the respective animal, and will have the same length.

# number of simulations per animal
nSimulations <- 5   
# perform simulations
simulations <- simulateMovement(trajectories = movementTrajectoriesCleaned, nSimulations = nSimulations, proj4 = proj4zone32)

Let’s visualize some simulations on an interactive map!

# plot example of simulated trajectory
exampleSim1 <- simulations[[2]][[1]][[1]]
exampleSim2 <- simulations[[2]][[2]][[1]]
# create spatial line
simExample1WGS84 <- createSpatialLine(locationdata = exampleSim1, proj4dataframe = proj4zone32, proj4line = proj4WGS84, lineId = animal)
simExample2WGS84 <- createSpatialLine(locationdata = exampleSim2, proj4dataframe = proj4zone32, proj4line = proj4WGS84, lineId = animal)
m <- leaflet() %>% 
  setView(lng=11.0367, lat=46.0331, zoom = 11) %>%
  # baselayers
  addTiles(group = 'Street map') %>%  
  addProviderTiles('Stamen.Terrain', group = 'Terrain') %>% 
  addProviderTiles('Esri.WorldImagery', group = 'Satellite') %>%
  # overlay groups
  addPolylines(data = simExample1WGS84, weight = 1, opacity = 0.8, color = "orange", group = 'Simulation 1') %>% 
  addPolylines(data = simExample2WGS84 , weight = 1, opacity = 0.8, color = "red", group = 'Simulation 2') %>% 
  addPolylines(data = trajectAlessandra, weight = 1, opacity = 0.8, color = "blue", group = 'Alessandra') %>% 
  addLayersControl(
    baseGroups = c('Street map', 'Terrain', 'Satellite'),
    overlayGroups = c('Alessandra', 'Simulation 1', 'Simulation 2'),
    options = layersControlOptions(collapsed = FALSE)
  )
m 

We extract the altitude and landuse values of simulated trajectories and store them in a matrix in order to compare them with the recorded GPS trajectories of our animals.

# extract altitude and landuse values of simulations 
simulationElevations <- extractRasterValues(DEM, simulations)
simulationLandUse <- extractRasterValues(landUse, simulations)
# combine recorded and simulated height data in one matrix
simElCombi <- combineSimTrack(heightDEM, animalNames, simulationElevations)
simLandCombi <- combineSimTrack(landuseTracks, animalNames, simulationLandUse)

Comparison recorded and simulated trajectories

Elevation

We can explore whether our roe deer have a preference for altitude through a visualization of movement on a DEM.

# variables
colors <- c("red","blue", "orange", "yellow", "purple")
names <- list(c("Agostino", "Allessandra", "Daiela", " Decimo", "Sandro"))
# select theme to plot with right colors
myTheme=rasterTheme(region=rev(brewer.pal('YlGn', n=9)))
# create plot of DEM with added trajectories
levelplot(DEM, main = "Altitude of recorded roe deer locations", margin = FALSE, par.settings=myTheme, xlim = c(647000, 670000), key = list(space = 'bottom', points = list(lty = 1, col = colors), text = names )) + 
    latticeExtra::layer({
  trajectAgostinoZone32 <- spTransform(trajectAgostino, CRS(proj4zone32))
  sp.lines(trajectAgostinoZone32, col = 'red')}) +
    latticeExtra::layer({
  trajectAllessandraZone32 <- spTransform(trajectAlessandra, CRS(proj4zone32))
  sp.lines(trajectAllessandraZone32, col = 'blue')}) +
    latticeExtra::layer({
  trajectDanielaZone32 <- spTransform(trajectDaniela, CRS(proj4zone32))
  sp.lines(trajectDanielaZone32, col = 'orange')}) +
    latticeExtra::layer({
  trajectDecimoZone32 <- spTransform(trajectDecimo, CRS(proj4zone32))
  sp.lines(trajectDecimoZone32, col = 'yellow')}) +
    latticeExtra::layer({
  trajectSandroZone32 <- spTransform(trajectSandro, CRS(proj4zone32))
  sp.lines(trajectSandroZone32, col = 'purple')})

Plots allow us to compare the recorded altitude of our animals with the simulations.

# plot the profiles of the simulated profiles and the trajectory
plotProfiles(simElCombi, animalNames)

The recorded trajectories show clearly less variation than the simulated trajectories. This result suggest that roe deer like to stay on a similar altitude for a longer timeperiod.

In these boxplots we can compare the distribution of the elevation of the recorded trajectories and the simulated trajectories to confirm our visual interpretation of the plotted profiles.

# variables for boxplot 
xlab <- "Trajectory"
ylab <- "Elevation (m)"
title <- "Comparison of recorded and simulated GPS trajectories of"
# visualize the difference between recorded and simulated data in boxplots
boxplotlist <- boxplotProfiles(simElCombi, animalNames, xlab, ylab, title, addName = TRUE) 
# arrage boxplots
grid.arrange(boxplotlist[[1]],boxplotlist[[2]],boxplotlist[[3]],boxplotlist[[4]],boxplotlist[[5]], layout_matrix = rbind(c(1,2), c(3,4),c(5,6)))

We can also compare the summary statistics (mean and standard deviation) of the recorded and simulated trajectories.

# calculate summary statistics for the recorded and simulated trajectories
summaryStats <- sumStats(simElCombi, animalNames)
# create empty vectors
xAxSim <- as.integer()
meansRec <- as.numeric()
meansSim <- as.numeric()
sdRec <- as.numeric()
sdSim <- as.numeric()
# store the summary statistics in a vector for plotting
for (animal in 1:length(animalNames)){
  meansRec <- c(meansRec, summaryStats[[animal]][1,1])
  meansSim <- c(meansSim, summaryStats[[animal]][1,2:ncol(summaryStats[[animal]])])
  sdRec <- c(sdRec, summaryStats[[animal]][3,1])
  sdSim <- c(sdSim, summaryStats[[animal]][3,2:ncol(summaryStats[[animal]])])
  
  # create an index vector for the x-axis
  for (col in 1:length(2:ncol(summaryStats[[animal]]))){
    xAxSim <- c(xAxSim, animal)
  }
}
# plot
par(mfrow=c(1,2), xpd = TRUE)
plot(xAxSim, meansSim, col = c('grey'), ylim = c(400,1600), main = 'Mean elevation', xlab = 'animal number', ylab = 'elevation (m)', pch = 19)
points(meansRec, col = 'red', pch = 19)

plot(xAxSim, sdSim, col = c('grey'), ylim = c(100,650), main = 'Standard deviation elevation', xlab = 'animal number', ylab = 'elevation (m)', pch = 19)
points(sdRec, col = 'red', pch = 19)

legend(-0.8,0, c("Recorded", "Simulated"), col = c("red", "grey"), pch = 19, bty = 'n')
par(mfrow=c(1,1), xpd = FALSE)

These various results all showa similar trend; roe deer stick to a relatively similar altitudes compared to the simulated trajectories. As well they seem to prefer higher altitude over lower altitude. However, this can be the result from human preference for lower altitudes…

Landuse

The next section analyses whether the roe deer have a preference for certain land use types. Please note that the following plot has double labels and a sub-optimal visualization. Sadly, we did not have enough time to figure these issues out. We hope you will forgive us these flaws!

# variables
colorsAnimals <- c("red","blue", "orange", "yellow", "purple")
colorsLandUse <- c('white','azure3', 'azure3', 'azure4','azure4','bisque4','bisque4','chartreuse', 'chartreuse2','chartreuse2','chartreuse3','darkolivegreen1', 'darkolivegreen1','darkolivegreen','darkolivegreen','darkorange2','darkorange2','darkgoldenrod3','darkcyan', 'darkslategray1')
names <- list(c("Agostino", "Allessandra", "Daiela", " Decimo", "Sandro"))
# copy landuse
landUseTest <- landUse 
# start for reclassification table
from <- 1
to <- as.integer()
start <- 1
clas <- as.integer()
# make reclassification table to remove double lables
for (label in 1:(length(landUseLegend$LABEL2)-1)){
  
  thisone <- landUseLegend$LABEL2[label]
  nextone <- landUseLegend$LABEL2[label+1]
  
  # check if label is same as next one
  if (!thisone == nextone){
    to <- c(to,label)
    from <- c(from, label+1)
  }
  
}
# add last number to get same length
to <- c(to,length(landUseLegend$LABEL2))
# create the reclassificatin table
rcl <- cbind(from, to, to)
# reclassify
landUseTest <- reclassify(landUseTest, rcl)
# ratify to get labels
landUseTest <- ratify(landUseTest)
rat <- levels(landUseTest)[[1]]
rat2 <- left_join(rat, landUseLegend, by = c("ID" = "GRID_CODE"))
rat2drop <- rat2[,c("ID","LABEL2")]
levels(landUseTest) <- rat2drop
# plot
levelplot(landUseTest, col.regions= colorsLandUse, opacity = 0.3, main = "Landuse of recorded roe deer locations", margin = FALSE, par.settings=myTheme, xlim = c(647000, 670000), key = list(space = 'bottom', points = list(lty = 1, col = colorsAnimals), text = names )) + 
    latticeExtra::layer({
  trajectAgostinoZone32 <- spTransform(trajectAgostino, CRS(proj4zone32))
  sp.lines(trajectAgostinoZone32, col = 'red')}) +
    latticeExtra::layer({
  trajectAllessandraZone32 <- spTransform(trajectAlessandra, CRS(proj4zone32))
  sp.lines(trajectAllessandraZone32, col = 'blue')}) +
    latticeExtra::layer({
  trajectDanielaZone32 <- spTransform(trajectDaniela, CRS(proj4zone32))
  sp.lines(trajectDanielaZone32, col = 'orange')}) +
    latticeExtra::layer({
  trajectDecimoZone32 <- spTransform(trajectDecimo, CRS(proj4zone32))
  sp.lines(trajectDecimoZone32, col = 'yellow')}) +
    latticeExtra::layer({
  trajectSandroZone32 <- spTransform(trajectSandro, CRS(proj4zone32))
  sp.lines(trajectSandroZone32, col = 'purple')})

The above map does not give us a very clear impression as to whether our animals prefer certain types of landuse. However, it does clearly show that the homerange of our animals mostly consists of forest and pasture.

Barplots may provide a clearer indication as to whether our animals have a preference for particular types of landuse:

# add legend to land use matrix for recorded and simulated data
simLandCombiDFLegend <- addLegend(simLandCombi, legend = landUseLegend, animalNames = animalNames)
# loop over the matrices for the different animals
for (animal in 1:length(animalNames)){
 
  # group based on ind (recorded / number of simulations)
  groupedLandUseAnimal <- as.data.frame(simLandCombiDFLegend[[animal]])
  groupedLandUseCount <- groupedLandUseAnimal %>% group_by(ind, LABEL2) %>% summarise(count = n())
  
  # create histogram
  assign(paste0("h",animal),ggplot(data=groupedLandUseCount, aes(x=LABEL2, y=count, fill=ind)) +
    geom_bar(stat="identity", position=position_dodge()) + ggtitle(paste("Land use of", animalNames[animal])) + xlab(NULL) + theme(axis.text.x = element_text(angle = 11, hjust = 1)))
}
grid.arrange(h1,h2, layout_matrix = rbind(c(1,1),c(2,2)))

grid.arrange(h3,h4, layout_matrix = rbind(c(1,1),c(2,2)))

grid.arrange(h5, layout_matrix = rbind(c(1,1),c(2,2)))

These plots confirm that our animals mostly reside in forest and pasture areas. Our simulated trajectories provide similar results to the recorded movement, which may explained by the fact that these landuse classes are prevalent in the area in which our animals live. The fact that this is their preferred home range does suggest that our roe deer love their strolls through forest and pasture!

Conclusion

This project has compared recorded and simulated roe deer movement in the area of Trento, Italy. It has analysed to what extent elevation and landuse can explain roe deer movement in the area. The results of the analysis suggest that the roe deer stick to a relatively small home range and stable altitude compared to the simulated locations. Roe deer had a preference for forest and pasture. However, the difference with simulated data is small, due to the high prevalence of these land use types in the home range of the roe deer.

# disconnect database
dbDisconnect(movementDB)

References

Urbano, F., & Cagnacci, F. (2014). Spatial database for GPS wildlife tracking data. Springer International Publishing.